Quantum factorization using three coupled harmonic oscillators 
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Using three coupled harmonic oscillators, we propose a quantum approach without using qubits 
to efficiently factorize an integer with only three steps. These three harmonic oscillators are coupled 
together via nonlinear interactions. This factoring method can be applied even if two harmonic 
oscillators are prepared in mixed states. As simple examples, we show how to factorize N — 15 and 
35 using this approach. For finding factors, the required time evolution of the system exponentially 
shortens for stronger nonlinear strengths and for higher-order nonlinearity. The effect of dissipation 
of the harmonic oscillators on the performance of this method is also discussed. Our results show 
that an alternative form of quantum computing can be exploited for factorization, without using 
qubits. 

PACS numbers: 03.67.Ac 
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Introduction. — Current cryptosy stems, such as RSA, 
depend on factorization of a large integer into a product 
of distinct prime numbers. In 1994, Peter Shor proposed 
[l| a quantum algorithm to factorize an integer in polyno- 
mial time. This algorithm gives an exponential speedup 
on the best known classical algorithm [2|, |3| . Implemen- 
tations of Shor's algorithm have been demonstrated in 
NMR [4[ and photonic systems [5|,|6[ to factorize N = 15. 
Recently, a new quantum algorithm [7] based on the adi- 
abatic theorem has been proposed to factorize integers, 
and it has been implemented (for N = 21) in an NMR ex- 
periment [7| . It is an interesting question to ask whether 
an alternative approach exists, to efficiently perform fac- 
torization. 

In this paper, we propose a new procedure to factor- 
ize integers by using three coupled harmonic oscillators, 
instead of qubits. In this method, an integer n is repre- 
sented as a number state \n) of a harmonic oscillator. The 
states of the first two harmonic oscillators are prepared 
in a number-state basis, which encode the trial factors 
of the number N to be factorized. The third oscillator 
is prepared in a coherent state. We will show that this 
approach without qubits can be applied even if the first 
two harmonic oscillators are in mixed states. In fact, 
the quantum coherence of the third harmonic oscillator 
is essential in this algorithm. This harmonic-oscillator 
quantum "computer" exhibits robustness against deco- 
herence. 

Now we outline the basic idea of this approach. We 
consider a model with specific nonlinear interactions be- 
tween three harmonic oscillators. In the phase-space rep- 
resentation, such nonlinear interactions make a coherent 
state rotate with a "phase angle" for a time t, where this 
coherent state and the number states of the two trial 
factors are in product states. This "phase angle" is pro- 
portional to the product of these two trial factors. The 
third oscillator thus acts as a "marker" for factors and 
non- factors. By quantum parallelism [8], the "phase an- 
gles" can be simultaneously computed for all trial factors. 



By performing a conditional measurement of a coherent 
state with a "phase angle" which is proportional to the 
product of two factors (i.e., the number N to be factor- 
ized), the resulting state of the oscillators 1 and 2 is the 
state of the factors. 

Three coupled harmonic oscillators. — We employ three 
coupled harmonic oscillators to perform quantum factor- 
ization. Its Hamiltonian is 
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H = ^y^cjjqtqj + h^^ 9k(ci\aia\a2) k ' a\a$, (1) 



i=i 



fc=i 



where ctj and ujj are the annihilation operator and the 
frequency of the j-th harmonic oscillator, respectively. 
The parameter g^ is the coupling strength of the non- 
linear interaction (a[aiala2) k a\a^ which is the product 
of the number operators of all three harmonic oscilla- 
tors, and the exponent k is a positive integer. Obviously, 
this number-conserving Hamiltonian H is exactly solv- 
able. The product of number states of the three harmonic 
oscillators |n,ra, I) = |n)i|ra)2|03 is an eigenstate of H 
with eigenvalue hyujin + uj2m + uj^I + X}/c=i( nm )^] • 

Quantum factoring. — Now we present a method to fac- 
torize an integer N. This approach involves only three 
steps: initialization, time evolution, and conditional mea- 
surement. Figure [T] shows the quantum circuit for this 
factoring method. 

1. Initialization. — The first two harmonic oscillators 
are prepared in either pure or mixed states, and the third 
harmonic oscillator is initialized as a coherent state, i.e., 



p(0)=pi(0)®p 2 (0)®p 3 (0), (2) 

Pi(0)®p2(0)= J2 P^ / |n,m>(n , ,m , |, (3) 



n,n' ,m,m' 



P3(0) = |a>3 3(a|, 



(4) 



where p™™ are the probabilities of the number states 
|n, ra)(n', m'\ of the oscillators 1 and 2, and n, n', m, m' = 
2, . . . , \N/2\. The trial factors are encoded to the state 
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FIG. 1. (Color online) Quantum factorization us- 

ing three harmonic oscillators. We prepare an input 
state p(0) = pi(0) ®p 2 (0) ®ps(0), where pi(0) p 2 (0) = 
E n , m Kr7T |n, m)(n / ,m / |, and p 3 (0) = |a) 3 3(a|. After ap- 
plying a unitary operator U(t), the reduced density matrix p r 
can be obtained by a conditional measurement of the coherent 
state |ajv(t))3, which rotates in phase space with a frequency 
Qzv. In the large- \a\ limit, the reduced density matrix p r 
becomes the state pf for the factors. 



of the harmonic oscillators 1 and 2, where each number 
state represents a trial factor. The eigenvalues nm of the 
state \n, m) are the product of the two trial factors n and 
m. 

2. Time evolution. — Let us now study the time 
evolution of the system starting with a product state 
\n,m)\a)3. The time-evolution operator U(t) = 
ex-p(-iHt) transforms the state \n^m)\a)s into 

U(t)\n,m)\a) 3 = e-^ n+ ^^|n,m)|c nm (t)) 3 , (5) 

where a nm (i) = exp(-iQ nm t)a, and Q nm = cu 3 + 
Efc=i^( nm ) fc i s t ne rotation frequency, in phase space, 
of the coherent state \a nm (t)). Here we have used 
ex.p(i f dalas)\a)s = |aexp(z#))3, where $ is a real num- 
ber [9[. The rotation frequency of the coherent state 
|<^nm)3 is increased by J2k=i 9k{nm) k , due to the non- 
linear interaction. If the product nm is equal to TV, 
the coherent state |a/v(£))3 rotates with a frequency 
Qn = ^3 + ^2k=i 9kN k in phase space. Otherwise, the 
coherent states for nm ^ N rotate with frequencies ^ nm , 
which are different to the frequency Q n. This means that 
the coherent states corresponding to factors' states and 
non- factors' states have different rotation frequencies in 
phase space. Thus, the state of the harmonic oscillator 
3 acts as a "marker" for the states of factors and non- 
factors. 

Now we apply the time-evolution operator U(t) to the 
initial state p(0). We then write the density matrix p(t) 
of the system as 

71,771,71' ,772/ 

(6) 

Where p%™ = e *[(n -n)wi + (m -m)w 2 ]t p n^m ? and n? m? n ' 

and m' denote all trial factors. By quantum parallelism 
[8], the functions a nm (t) are simultaneously computed 
for different values of n and m. Obviously, this state 



becomes non-separable [10| between the three harmonic 
oscillators. The quantum entanglement [ll| between the 
harmonic oscillators is produced due to the nonlinear in- 
teractions of the three harmonic oscillators. 

3. Conditional measurement. — We then perform a 
measurement conditional on the coherent state \aj\r (£))$, 
where this coherent state and the factor's state are in 
product states at the time t. In this way, we can 
abandon the state of non-factors if the orthogonality 
of coherent states is approximately valid. Here we de- 
note the products rs and rV to be equal to N (i.e., 
rs,r's' = N). The resulting reduced density matrix 
p r = Tis{\cYN(t)}3 s(aN{t)\p(i)}/A can be obtained with 
probability J2 r s Prs as 



p r (t) 
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,|ra,ra)(ra',ra'|, (7) 



71,771,71' ,777' 



where A = Tr{|aAr(£))3 s(ajsf(t)\p{t)} is a normalization 
constant as ]C n ,m^nml e nm| 2 : , and the coefficient e nm = 
(a N (t)\a nm (t)} = exp{-\a\ 2 [l-exp(in N t-in nm t)}} 0> 
is the overlap between the two coherent states \ajsr(i)) 
and \a nm (t)). 

In the limit of large \a\, the orthogonality of coherent 
states holds approximately, such that the overlap e nrn 
tends to zero. Thus, we can obtain the state of the factors 



lim p r = p f , 

\a\— t-oo 

where pf is the density matrix of the factors 



(8) 



(9) 



r,s,r' ,s' 



and Af = ^2 r sPrs i s a normalization constant. Indeed, 
the reduced density matrix p r is well approximated by the 
state of factors p/, if \a\ is sufficiently large. Finally, we 
can determine the factors of the number N by measuring 
the state of the harmonic oscillators 1 and 2. 

We also remark that the overlap between the two co- 
herent states e nrn is a function of time t. To achieve the 
best approximation of the state of the two factors, the 
overlap between the coherent state \ajsr(t)) and the other 
coherent states in Eq. ([?)) can be minimized by perform- 
ing the conditional measurement at an appropriate time 
t*. This time t* depends on the initial probability distri- 
bution and the non-linear interaction strength. 

Example: Initial thermal states. — We consider the ini- 
tial state of the first two harmonic oscillators in ther- 
mal states and the third harmonic oscillator is a coherent 
state, i.e., 

P(°) = y^ J PinP2m\n,m)(n,m\ <g) |a) 33 (a|, (10) 



where p jn = [1 - e -^ 3 ^ e -pn 3 huj 3 ig the probability of 
the thermal states of the j-th harmonic oscillator [9[, 
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FIG. 2. (Color online) The fidelity F(t/u 3 ) between the 
reduced density matrix p r (t/ 003) and the factor's states pf for 
N — 15, is plotted as a function of the dimensionless time 
£/cj3, for the parameters cji = 1.5cj3, 0J2 — 2cj3, T — 3huo3/kB 
and \a\ — 5. The mean excitation numbers fij are 1.542, 
1.055 and 2.528, for j = 1,2,3 and k = 1. The probability 
of obtaining the coherent state \QLN=i5(t/u)3)) is 3.65 x 10 -3 . 
Black, blue and red lines correspond to the results for g — CJ3, 
0.9a;3, and 0.8cj3- For different strengths g, the maxima of 
the respective fidelities are indicated by the arrows. 



/3 = l//c#T, ks is the Boltzmann's constant, T is the 
temperature, and j = 1,2. The mean excitation number 
fij is [ e huj j/k B T — 1]- 1 for the bosonic mode cij [9]. Here 
we require that the probabilities of the number states 
|n, m) are much smaller than that of the state of factors, 
where n and m are larger than N. 

As a simple illustrative example, we now factor N = 
15 = 3 x 5 using this quantum approach. For simplicity, 
we consider the Hamiltonian as 



Hk = h2_.UjCLjaj + hg(a\aiala2) k alas 1 



(11) 



where g is the nonlinear strength and k is a positive in- 
teger. This allows us to clearly study the role of nonlin- 
earity in this method. 

To evaluate the performance of this approach, we in- 
vestigate the fidelity between the reduced density matrix 
p r (t) and the state of the factors pf at the time t. The 



fidelity is defined as [12 1 



Fw={Tv[„yv«) ( .y 2 ]" 2 



r 



(12) 



The fidelity F(t) = ^2 r ^ s PirP2s/ ^Z n ,mP^riP2m\^nm\ 2 can 
be obtained for the initial thermal states in Eq. (fT0]h 
where r, s = 3, 5. We now take the exponent k in Eq. 
(ITT]) to be one, i.e., k = 1. In Fig. [2j the fidelity F(t/uj 3 ) 
of the state of the factors is plotted as a function of the 
dimensionless time t for the magnitude \a\ = 5 and differ- 
ent values of the dimensionless nonlinear strengths g/uJs. 
The fidelity F begins to increase as a function of the time 
and reaches a maximum at t* /us « 0.335, for g = 003. 
Then, the fidelity sharply drops and fluctuates with time 
t* /0J3. Fig. [2] shows that the fidelities exhibit similar 
patterns versus t/us for different nonlinear strengths g. 



FIG. 3. (Color online) The logarithmic dimensionless time 
]n(t/uj3) is plotted as a function of the nonlinear strength 
g/003 for \a\ =6 and the same parameters as in Fig. O The 
time t/uj3 is taken when the fidelity F for N — 15 first arrives 
at 0.9. Black, red, blue, and green dots correspond to the 
results for k — 1, 2, 3, and 4. 



Note that, for the smaller strength g, the fidelity F takes 
longer time to reach the maximum (F c^ 1). 

We now proceed to study the performance of the al- 
gorithm with higher-order k nonlinearity and stronger 
nonlinear interaction strength g. In Fig. [3j we plot the 
logarithmic dimensionless time In (t/us) as a function of 
g for I a I = 6 and k — 1. The time i/uj$ is taken when 
the fidelity F first arrives at 0.9. We can see that the 
time In (i/us) steadily decreases as g/uus increases. In 
the same figure, we also show results for the higher val- 
ues of k (=2,3 and 4). The logarithmic time decreases 
by 3 as k increases by 1. Thus, a higher-order nonlin- 
earity can exponentially shorten the evolution time to 
achieve the high fidelity F* > 0.9. This result shows 
that both the high- order k nonlinearity and strong non- 
linear strength g are very useful to increase the efficiency 
for finding factors. 

In Fig. [4j we plot the fidelity F(t* /ujs) versus \a\, for 
t* /cjOs = 0.335 and g\ = U3. The fidelity F increases 
as the value of \a\ increases. The fidelity F becomes 
saturated and reaches near unity when \a\ ~ 7. Note 
that the fidelity cannot be exactly equal to one because 
the probabilities of the initial number states |n,m) are 
not exactly zero where n and m are larger than 15. 

Effect of dissipation. — We now consider the system 
weakly coupled to a thermal environment. We also as- 
sume that the Markovian approximation is valid [13j |. 
The master equation of the density matrix then reads 



-i[H,p] + JT t [£ j (p) + £i(p)]> 



(13) 



j=i 



where Cj(p) = Jj(rij + l)(2ajpo)- 



CLnalp 



(LjCLjP 



pa]a,j)/2, 



paja}j)/2, and jj is the 



Cj(p) = -i j n j (2a)pa j - <^ jr r ^ j; 
damping rate of the j-th harmonic oscillator. 

Now we consider the dissipative dynamics of the sys- 
tem for the initial state in Eq. ([TO]) . Since the states of 
the harmonic oscillators 1 and 2 are prepared in thermal 
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FIG. 4. (Color online) The fidelity F(t*/u> 3 ) for AT = 15, 
plotted as a function of \a\ at the time £*/cj3, for the nonlinear 
strength g = 003 and the same parameters as in Fig. [21 The 
fidelities vary with the dissipation rate 73 as shown. 



states which are the steady states of harmonic oscilla- 
tors in the thermal environment, we can write the total 
density matrix p(t) = ^ n ^ m PinP2 m \n,m)(n,m\(3p% rn (t), 
where ps m (t) is the density matrix of the harmonic oscil- 
lator 3 corresponding to the number states |n, ra)(n, m\. 
Only the harmonic oscillator 3 is subject to the dissipa- 
tion. This allows us to write the master equation of the 
density matrix p^ m as 



■nm 
P3 



-ift. 



nm [a 3 <23, p™] 



■CM 



nm ) + £M m )- (14) 



This master equation is exactly solvable [14[. Therefore, 
the time evolution of the total density matrix p{t) can be 
solved by summing over all possible solutions Ps m (t) in 
Eq. (fl4|) corresponding to the probabilities p\ n P2m- 



strength g. However, we can see that the effect of dissi- 
pation on the performance of this factoring algorithm is 
negligible when \a\ is about 8. The high fidelity can be 
attained in the presence of dissipation if the value of \a\ 
is sufficiently large. 

In addition, this approach can be used for factoring 
larger integers. We show another example for factoring 
N = 35 = 5 x 7. In Fig. El we plot the fidelity F(t/uj 3 ) of 
the factor's states for N = 35 versus the time t/uo^. High 
fidelities can be obtained by increasing the magnitude 
\a\ but the probability of obtaining the coherent state 
\<^N=35(t/oJs)} becomes very low. This probability can 
indeed be increased by appropriately choosing an initial 
thermal state. 

Conclusion. — We have proposed a non-qubit-based ap- 
proach to factor an integer using three coupled harmonic 
oscillators. This state of factors of an integer can be ob- 
tained by only three steps. This method can also work if 
the harmonic oscillators 1 and 2 are in mixed states. The 
efficiency for finding factors significantly increases if a 
higher-order nonlinearity and stronger nonlinear strength 
are used. The effect of dissipation on the performance of 
this factoring approach was also studied. 
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FIG. 5. (Color online) The fidelity F{t/m) for N = 35, 
plotted versus the time t/uj% , for g = 003 and the same param- 
eters as in Fig. [2] The probability of obtaining the coherent 
state \aN=35(t/u)3)) is 3.54 x 10 -4 . Black, blue, and red lines 
correspond to the results for \a\ =6, 8, and 10. 

In Fig. |U we plot the fidelity of the states for the fac- 
tors versus \a\ subject to the dissipations 73 = (black 
circle dotted line), 73 = 0.5co>3 (blue square-dotted line) 
and 73 = CJ3 (red star-dotted line). As shown in Fig. [4j 
the fidelity F(t*/us) decreases for small \a\ if the dissi- 
pation rate 73 is comparable to the nonlinear interaction 
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